home *** CD-ROM | disk | FTP | other *** search
/ Precision Software Appli…tions Silver Collection 4 / Precision Software Applications Silver Collection Volume 4 (1993).iso / stats / chadyn.exe / YBASINS.C < prev    next >
Text File  |  1988-12-15  |  34KB  |  1,133 lines
  1.  
  2. /***************************  YBasins.c  *************************************/
  3. /********************* (C) 1986,7,8 by JAMES A. YORKE ************************/
  4. /********** ROUTINES in file:
  5. SetNum_tests() sets num_tests which is the number of points to be checked by 
  6.         checkTrajectory; we will not check points if a coordinate (the
  7.         zeroth coordinate) is -9999; we want num_tests-1 to be the 
  8.         last storage vector for which that coordinate is not -9999  
  9. setframe(lower,upper,pedges) given two numbers, lower and upper, a new 
  10.         interval (larger than the original if diameters[ScrnSec] > 0 ) 
  11.         is found and the ends are stored in pedges[0] and pedges[1] 
  12.         where pedges[0] < pedges[1] 
  13.     double lower,upper,*pedges;
  14. BasinBoundary()   initial points are plotted when storage vector 2 is not 
  15.         set and the trajectory goes to infinity , that is when
  16.             (tau == 1 && y[eqn0 + eqn1 + zeroth] == -9999)
  17.         OR when the trajecory approaches storage vector 21  
  18. int checkTrajectory(EnterLevel,pVec)
  19.         This tests where the trajectory-that-starts-at-y is 
  20.         initially and where it is after 
  21.         each call of (*ChkTrajIteratee)(); checking the trajectory 
  22.         for a maximum of "MaxChecks" times. The routine leaves 'time' 
  23.         equal to the exit time (more or less). The routine:
  24.  
  25.         returns 0 if it has gone thru the full set of MaxChecks checks
  26.             without getting near any of the specified attractors 
  27.             and without leaving the box 
  28.         returns 1 if it is ever outside the box frame[4];
  29.         returns k (2 <= k <= 7) if the trajectory of the point 
  30.             approaches close to point yk[]  /
  31.     int EnterLevel;
  32.     double *pVec;
  33. int check_coord(pleft,Y)
  34.          pleft is a pointer to the value of the coordinate for 
  35.         the left side of the box. We assume pleft+1 points to the 
  36.         value of the right side. this function checks the two sides of 
  37.         the box for a given number Y. This function returns 
  38.                 IN = 1 if point y is in interval,  
  39.                 OUT = 0 if outside, and 
  40.                 NOTSET = -9999 if neither box coord is defined. /
  41.     double *pleft,Y;
  42. ************************* STRADDLE ORBIT ROUTINES ************************
  43. initSadStraddle()/* first we find out if ya and yb are set 
  44.         and if so, what checkTrajectory returns for ya and yb; 
  45.         the returned values are retA and retB; 
  46.         if retA = retB, or if either equals 0, quit by decreasing 
  47.         level; /
  48.  
  49. SetBasin() sets the vector basin[8] so that each coordinate is either 
  50.         1 or 2 thereby indicating which attactor each yx[] point 
  51.         is in; hence when getting a number that is returned from 
  52.         checkTrajectory, the program can tell which basin it is in /
  53. initBasStraddle() first we find out if ya and yb are set 
  54.         and if some numbered storage vector >= 2 is set
  55.         and if so, what checkTrajectory returns for ya and yb; 
  56.         the returned values are retA and retB; 
  57.         if retA = retB, or if either equals 0, quit by decreasing 
  58.         level; /
  59. DoAB_Exit() checks where ya and yb are going and complains if they go 
  60.         nowhere /
  61. int AreAB_set() if ya and yb are set 1 is returned; otherwise 0 /
  62. initChTraj()
  63. BasinStraddle()
  64.      this routine takes the pair of vectors ya,yb and refines the pair
  65.         until they are "IsClose"; then it takes a single iterate of
  66.         the refined ya and yb so that the pair advances one iterate;
  67.         the method of refinement depends on BST_switch: 
  68.         BST_switch = 1 means a random point between 
  69.         ya and yb is chosen; otherwise the midpoint is chosen /
  70. AccessBS() AccessBasinStraddle: this routine 
  71.         takes the pair of vectors ya,yb and refines the pair
  72.         until they are "IsClose"; then it takes a single iterate of
  73.         the refined ya and yb so that the pair advances one iterate;
  74.         it is refined by taking lots of points between ya and yb and 
  75.         as soon as any of them are found to be in basinB, move yb to 
  76.         that point. If none of them are found to be in basinB, choose 
  77.         a point near the midpoint that is known to be a basinA point 
  78.         and move ya to it  /
  79. SaddleStraddle() we search for a trajectory that stays away from the basins;
  80.          this routine takes the pair of vectors ya,yb and refines the pair
  81.         until they are "IsClose"; then it takes a single iterate of
  82.         the refined ya and yb so that the pair advances one iterate /
  83. int refineSadStrad()
  84.          this routine tries to find a refinement, i.e. a new pair
  85.         ya, yb, (one could be unchanged) closer together than the 
  86.         starting pair; if it is successful, the pair is changed and
  87.         1 is returned; 
  88.         it returns 0 if it finds an orbit that never leaves the domain;
  89.         it returns -1 if it cannot find a refinement; notice that all 
  90.         the values of fraction are diadic and so are computed exactly
  91. OneIterate(Y)   iterates the point at Y[] using y /
  92.     double *Y;
  93. set_y_fraction(fraction)
  94.     double fraction;
  95. GiveUp()
  96. triple()
  97. ******************************************************************************/
  98.  
  99. /* "never exits" case must be reworked to exit smoothly from program */
  100. #include "yinclud.h"
  101.  
  102. #define OUT 0
  103. #define IN  1
  104. #define NOTSET -9999
  105.  
  106.  
  107. double frame[4];
  108. int num_tests;
  109. static int isYkSet[9];
  110. static int time;
  111. static int retA,retB;
  112. static int basin[8];
  113. static int basinA,basinB;
  114. static int intrptFlag = NO;
  115. static double totalRefinementCalls,totalTime;
  116. int modNum = 14;
  117.  
  118.  
  119. SetNum_tests()/* sets num_tests which is the number of points to be checked by 
  120.         checkTrajectory; we will not check points if a coordinate (the
  121.         zeroth coordinate) is -9999; we want num_tests-1 to be the 
  122.         last storage vector for which that coordinate is not -9999  */
  123. {
  124.     int k;    
  125.     for (num_tests = NUM_Y;
  126.          y[eqn0 +(num_tests - 2)*eqn1 + zeroth] == -9999 && num_tests >= 2;
  127.          num_tests--)
  128.         ;
  129.     isYkSet[8] = NO;
  130.     for (k = 2; k <= 7; k++)
  131.         if (y[eqn0 + eqn1*(k-1) + zeroth] != -9999. )
  132.             isYkSet[k] = YES;
  133.         else
  134.             isYkSet[k] = NO;
  135. }
  136.  
  137.  
  138. setframe(lower,upper,pedges)/* given two numbers, lower and upper, a new 
  139.         interval (larger than the original if diameters[ScrnSec] > 0 ) 
  140.         is found and the ends are stored in pedges[0] and pedges[1] 
  141.         where pedges[0] < pedges[1] */
  142.     double lower,upper,*pedges;
  143. {
  144.     double temp;
  145.  
  146.     pedges[0] = lower - diameters[ScrnSec]*(upper - lower);
  147.     pedges[1] = upper + diameters[ScrnSec]*(upper - lower);
  148.  
  149.     if (lower > upper) /* switch so pedges[0] will be less than pedges[1]*/
  150.     {
  151.         temp = pedges[0];
  152.         pedges[0] = pedges[1];
  153.         pedges[1] = temp;
  154.     }
  155. }
  156. extern long initialdot;/* see Y.C */
  157. BasinBoundary()/* initial points are plotted when storage vector 2 is not 
  158.         set and the trajectory goes to infinity , that is when
  159.             (tau == 1 && y[eqn0 + eqn1 + zeroth] == -9999)
  160.         OR when the trajecory approaches storage vector 21  */
  161. {
  162.     long jy;
  163.     int tau, level0=level, oldColor = color;
  164.     int cheqFreq();
  165.     double x0 ,x1, y0, y1, alittle,timeofday();
  166.  
  167.     X_step = ( X_upper - X_lower)/xpixels;
  168.     Y_step = ( Y_upper - Y_lower)/ypixels;
  169.  
  170.     level++;
  171.     if (level == PROCESS)
  172.         boot_crt();/* sets boot to 0, and if boot is 0, it clears
  173.                 p[][] and the screen. It sets the scale */
  174.     SetBasin();/* sets basin[] */
  175.     initChTraj();
  176.  
  177.     if (VertLine > 0 && VertLine < xpixels && refreshFlag == YES)
  178.         ;
  179.     else
  180.         VertLine = 0;
  181.     refreshFlag = NO;
  182.  
  183.     if (TDFreq > 0)
  184.         PRINT 
  185.      "WARNING: picture is stored to disk every %ld vertical line(s) of dots.\n"
  186.                 ,TDFreq);
  187.  
  188.     time0 = timeofday();
  189.     dot = dotAtTime0 = ((long) VertLine)*1000;
  190.     initialdot = ((long) VertLine)*1000;
  191.  
  192.  
  193.     initTime();
  194.        for (;  (level == level0+1)  &&  VertLine < xpixels; VertLine++)
  195.        {    
  196.         if(timeFlag == YES || TDFreq > 0)
  197.             if(checkFreq() == YES) 
  198.             {    /* means terminate; checkFreq also handles
  199.                    occasional transmissions of picture to disk 
  200.                 */
  201.                 PRINT "process termination: OUT OF TIME \n");
  202.                 level = level0 - 1;
  203.                 continue;
  204.             }
  205.  
  206.         alittle = .000001;
  207.         x0 = X_lower + X_step*(VertLine + alittle);
  208.         x1 = X_lower + X_step*(VertLine + 1 + alittle);
  209.         dot = VertLine*1000 ;/* 3 least significant digits = height */
  210.         for (jy = 0; jy < ypixels && level == level0 + 1; jy++)
  211.         {
  212.             dot++;
  213.             y0 = Y_upper - Y_step*(jy + alittle);
  214.             y1 = Y_upper - Y_step*(jy + 1 + alittle);
  215.             if (dim > 2) setequal(0,1,eqn0);
  216.             y[X_coord] = x0;
  217.             y[Y_coord] = y0;
  218.  
  219.             tau = checkTrajectory(level,y);
  220.  
  221.                if ( basin[tau] == 1)
  222.                {
  223.                 color = (time % modNum)+1;
  224.                 _setcolor(color);
  225.                 block_plot(x0,x1,y0,y1);
  226.                 /* in high resolution modes it is 
  227.                 possible that for some vertical lines of pixels
  228.                 x0 = x1 within the resolution of the screen so
  229.                 nothing will be plotted on the screen; it will
  230.                 however be plotting in p[][] */
  231.                }    
  232. #ifndef MAINFRAME
  233.             if (ChkKeyStatus() != 0 || cycle > 0) 
  234.                 Interrupt();
  235. #endif
  236.         }
  237.        }
  238.     if (VertLine != xpixels)
  239.         VertLine--;/* if we terminate in the middle of a picture, and 
  240.             we begin again later, we want to begin at the 
  241.             beginning of a line */
  242.     else
  243.     {
  244.         VertLine = 0;/* no harm to reset it if the entire picture 
  245.             has been redrawn */
  246.     }
  247.     if ((onprint>=1) && (level == level0+1))
  248.             pprint(onprint);/* this sends picture to B: disk if disk = 1 */
  249. #ifndef MAINFRAME
  250.     if (level == level0+1)
  251.     {
  252.         cycle = 1;
  253.         Interrupt();
  254.     }
  255. #endif
  256.     level = level0;
  257.     setequal(0,1,eqn0);/* reinitializes y[] = y1 */
  258.     if (level < PROCESS)
  259.     {
  260.         scr_clr();/* clears screen but not pic[]; it is a function 
  261.                             in desmets pcio.a */
  262.         MainMenu();
  263.     }
  264.     color = oldColor;
  265. }
  266.  
  267. whenAndWhere()/* tells where the trajectory through y1 goes and how long 
  268.         it takes to get there; for interrupt W and command W */
  269. {
  270.     double yTemp[EQUATIONS];
  271.     int ret,checkTrajectory();
  272.     int oldNumLyap = num_lyap;
  273.  
  274.     num_lyap = 0;
  275.     store(yTemp,y,lyapzero);/* this routine operates on y[] */
  276.  
  277.     SetBasin();/* sets basin[] */
  278.     initChTraj();
  279.  
  280.     ret = checkTrajectory(level,y+eqn0);/* test y1 */
  281.     store(y,yTemp,lyapzero);/* this routine operates on y[] */
  282.     num_lyap = oldNumLyap;
  283.  
  284.     if (ret == 0 && printer > 0)
  285.         PRINT 
  286. "Trajectory through small cross iterated MC=%d times, end point stored in ye\n"
  287.         ,MaxChecks);
  288.     if (ret == 1 && printer > 0)
  289.         PRINT 
  290. "Trajectory through small cross diverges toward infinity in %d iterates.   \n"
  291.         ,time);
  292.  
  293.     if (ret >= 2 && ret <= 7 && printer > 0)
  294.         PRINT 
  295. "Trajectory through small cross came within RA=%3.3lf of y%d in %d iterates\n"
  296.         ,ra[ret],ret,time);
  297. }
  298.  
  299. int checkTrajectory(EnterLevel,pVec)
  300.         /* This tests where the trajectory-that-starts-at-y is 
  301.         initially and where it is after 
  302.         each call of (*ChkTrajIteratee)(); checking the trajectory 
  303.         for a maximum of "MaxChecks" times. The routine leaves 'time' 
  304.         equal to the exit time (more or less). The routine:
  305.  
  306.         returns 0 if it has gone thru the full set of MaxChecks checks
  307.             without getting near any of the specified attractors 
  308.             and without leaving the box 
  309.         returns 1 if it is ever outside the box frame[4];
  310.         returns k (2 <= k <= 7) if the trajectory of the point 
  311.             approaches within RA of the point yk[]  */
  312.     int EnterLevel;
  313.     double *pVec;
  314. {
  315.     int k;
  316.     double distance(),*p1,*p2;
  317.  
  318.  
  319.        if (pVec != y)
  320.         store(y,pVec,lyapzero);/* this routine operates on y[] */
  321.  
  322.     p1 = y+zeroth;
  323.     p2 = y + eqn0 + zeroth;
  324.     for (time = 0; (time < MaxChecks) && (level == EnterLevel); time++)
  325.     {
  326.  
  327.            /* We now handle the case where the point is "diameters" outside
  328.               the screen region */
  329.         if( check_coord(frame,  y[X_coord]) == OUT || 
  330.             check_coord(frame+2,y[Y_coord]) == OUT )   
  331.         {
  332.             return (1);
  333.         }
  334.  
  335.         for (k = 2; k < num_tests; k++)
  336.             if (isYkSet[k] == YES )/* has yk[] been set? */
  337.             {
  338.                   if (distance(p1, p2 + (k-1)*eqn1, vec_dim) 
  339.                             < ra[k])
  340.                {
  341.                 return(k);/* Y is close to kth test point */
  342.                }
  343.             }
  344.         (*ChkTrajIteratee)();
  345.        }
  346.     if (time == MaxChecks && MaxChecks >= 100 && level > 1)
  347.     {
  348.         store(ye,y,lyapzero);
  349.         PRINT 
  350. "A trajectory has been iterated MC times; final point is stored in ye[]   \n");
  351.         PRINT 
  352. "Either set some y2...y7 to ye    or use command MC to set MC to be about 20\n"
  353.         );
  354.         
  355.     }
  356.        return (0);/* one way to get here is for Interrupt() to decrease
  357.         level via interrupt ; another way is for the trajectory to be 
  358.         iterated the full MaxCheck iterates without meeting any of
  359.         the tests; */
  360. }
  361.  
  362.  
  363. int check_coord(pleft,Y)
  364.         /* pleft is a pointer to the value of the coordinate for 
  365.         the left side of the box. We assume pleft+1 points to the 
  366.         value of the right side. this function checks the two sides of 
  367.         the box for a given number Y. This function returns 
  368.                 IN = 1 if point y is in interval,  
  369.                 OUT = 0 if outside, and 
  370.                 NOTSET = -9999 if neither box coord is defined. */
  371.     double *pleft,Y;
  372. {
  373.     int ret;
  374.  
  375.     if((*pleft == NOTSET) && (*(pleft+1) == NOTSET)) 
  376.         ret = NOTSET;
  377.     else 
  378.     {
  379.       ret = ((*pleft != NOTSET)&&(*pleft > Y) ? OUT : IN );
  380.       if (ret == IN)
  381.           ret = ((*(pleft+1) != NOTSET)&&(*(pleft+1)) < Y ? OUT : IN);
  382.     }
  383.     return (ret);
  384. }
  385.  
  386. /************************* STRADDLE ORBIT ROUTINES ************************/
  387.  
  388. initSadStraddle()/* first we find out if ya and yb are set 
  389.         and if so, what checkTrajectory returns for ya and yb; 
  390.         the returned values are retA and retB; 
  391.         if retA = retB, or if either equals 0, quit by decreasing 
  392.         level; */
  393. {
  394.     totalRefinementCalls = totalTime = 0;
  395.  
  396.     if ((AreAB_set()) == 0) /* 0 means>= one is not set;
  397.             if 0, then level is downgraded */
  398.         return;
  399.     initChTraj();
  400.     
  401.     DoAB_Exit();/* If checkTrajectory returns 0, then level is downgraded*/
  402. }
  403.  
  404. SetBasin()/* sets the vector basin[8] so that each coordinate is either 
  405.         1 or 2 thereby indicating which attactor each yx[] point 
  406.         is in; hence when getting a number that is returned from 
  407.         checkTrajectory, the program can tell which basin it is in */
  408. {
  409.     /* the number ret that checkTrajectory returns is:
  410.         0 if the trajectory never exits,
  411.         1 if it goes to infinity, i.e., leaves the box;
  412.         k if it enters the ball centered at yk[];
  413.  
  414.         now we say basin is #1 if ret = 2,3, or 4 and
  415.                  is #2 if ret = 5,6, or 7 or 0;
  416.         ret = 1 (infinity case: if it is ever outside the box 
  417.         frame[4]) is basin #1 IF y2[0] is default = -9999.;
  418.         otherwise ret = 1 is basin #2  */
  419.          
  420.  
  421.     basin[2] = basin[3] = basin[4] = 1;
  422.     basin[0] = basin[5] = basin[6] = basin[7] = 2; 
  423.     if (y[eqn0+eqn1] == -9999.)
  424.         basin[1] = 1;
  425.     else
  426.         basin[1] = 2;
  427. /*
  428. PRINT 
  429. "basin[0,7] = %d    %d    %d %d %d    %d %d %d ]\n"
  430. ,basin[0],basin[1],basin[2],basin[3],basin[4],basin[5],basin[6],basin[7]);
  431. */
  432. }
  433.  
  434. initBasStraddle()/* first we find out if ya and yb are set 
  435.         and if some numbered storage vector >= 2 is set
  436.         and if so, what checkTrajectory returns for ya and yb; 
  437.         the returned values are retA and retB; 
  438.         if retA = retB, or if either equals 0, quit by decreasing 
  439.         level; */
  440. {
  441.     int initLevel = level;
  442.  
  443.     SetBasin();
  444.     if (AreAB_set() == 0) /* 0 means>= one is not set */
  445.         return;
  446.     initChTraj();
  447.     DoAB_Exit();/* defines retA and retB */
  448.     basinA = basin[retA];
  449.     basinB = basin[retB];
  450.  
  451.     initChTraj();
  452.     if (num_tests <= 2)/* set in initChTraj() */
  453.     {
  454.         PRINT
  455.               "\n   @#$%^&*+#    THERE IS A PROBLEM   *&%$#*|@+ \n\n");
  456.         PRINT
  457. "You must first tell the program where at least one attractor is. You must set"
  458.         );
  459.         PRINT
  460.            "\n            at least one of the vectors 2 thru %d.\n",NUM_Y);
  461.         level = initLevel-1;
  462.         return;
  463.     }
  464.     
  465.     if (basin[retA] == basin[retB])/* are they in the same basin? */
  466.     {
  467. PRINT 
  468. "retA=5d; retB=%d\n",retA,retB);
  469. PRINT 
  470. "basin[retA]=5d; basin[retB]=%d\n",basin[retA],basin[retB]);
  471.  
  472.         PRINT
  473. "Verrry baaad: stored points a & b are in the same basin %d\n"
  474.             ,basin[retA]);
  475.         PRINT
  476. "This occurs sometimes when the radius RA is too big, giving disks \n"
  477. );
  478.         PRINT 
  479. "        that actually lie in two basins\n");
  480.         level = initLevel-1;
  481.     }
  482.     store(y,yb,lyapzero);/* set y to equal yb so that the first point 
  483.             plotted will not be a bad point; do not change 
  484.             high coordinates of ya needed for exponents */
  485. }
  486.  
  487.     
  488. DoAB_Exit()/* checks where ya and yb are going and complains if they go 
  489.         nowhere */
  490. {
  491.     int initLevel = level;
  492.  
  493.     storeNumLyap = num_lyap;
  494.     num_lyap = 0;/* do not compute exponents during this phase */
  495.  
  496.     retA = checkTrajectory(level,ya);
  497.     if (printer >= 2)
  498.         PRINT
  499.            "ya returns attractor# %d in time %d                 \n"
  500.             ,retA,time);
  501.     retB = checkTrajectory(level,yb);
  502.     if (printer >= 2)  
  503.         PRINT
  504.         "yb returns attractor# %d in time %d                      \n\n"
  505.             ,retB,time);
  506.     num_lyap = storeNumLyap;
  507.  
  508.     if (retA == 0) 
  509.     {PRINT
  510.     "stored point a returns zero: its trajectory goes to no attractor \n");
  511.         level = initLevel - 1;
  512.         return;
  513.     }
  514.     if (retB == 0) 
  515.     {PRINT
  516.     "stored point b returns zero: its trajectory goes to no attractor \n");
  517.         level = initLevel-1;
  518.         return;
  519.     }
  520. }
  521.  
  522.  
  523. int AreAB_set()/* if ya and yb are set 1 is returned; otherwise 0 */
  524. {
  525.     if (ya[zeroth] == -9999)
  526.     {   PRINT
  527.         "stored point a is still at the default value; cannot proceed \n");
  528.         level--;
  529.         return(0);
  530.     }
  531.     if (yb[zeroth] == -9999)
  532.     {   PRINT
  533.         "stored point b is still at the default value; cannot proceed \n");
  534.         level--;
  535.         return(0);
  536.     }
  537.     return (1);
  538. }
  539.  
  540. initChTraj()
  541. {
  542.     int k;
  543.        /* checkTrajectory checks to see if the trajectory has left the box with
  544.         edges in frame[] so we have so set frame[] */
  545.     setframe(X_Lo[ScrnSec],X_Up[ScrnSec],frame);
  546.     setframe(Y_Lo[ScrnSec],Y_Up[ScrnSec],frame+2);/* the box with 
  547.         coordinates 
  548.         frame[0] < x < frame[1]; frame[2] < y < frame[3] 
  549.         is larger than the screen box by a factor 1 + 2* diameters */
  550.     SetNum_tests();/* num_tests is the number of points to be checked by 
  551.       checkTrajectory; we will not check points if a coordinate is -9999 */
  552.  
  553.     for (k=2; k<9; k++)
  554.         ra[k] = rad_attr;
  555.     if (ra2 != -9999.)/* ra2 and ra5 are set using commands RA2 & RA5 
  556.                 in YCOMB.C; see also YPICKMAP.C*/
  557.         ra[2] = ra2;
  558.     if (ra5 != -9999.)
  559.         ra[5] = ra5;
  560. }
  561.  
  562.  
  563. BasinStraddle()
  564.     /* this routine takes the pair of vectors ya,yb and refines the pair
  565.         until they are "IsClose"; then it takes a single iterate of
  566.         the refined ya and yb so that the pair advances one iterate;
  567.         the method of refinement depends on BST_switch: 
  568.         BST_switch = 1 means a random point between 
  569.         ya and yb is chosen; otherwise the midpoint is chosen */
  570. {
  571.     int random997(); /* in YTOOLS.C Returns an integer between 0 and 997 */
  572.     int ret;
  573.     int initLevel = level;
  574.     double dist,distance(),fraction;
  575.  
  576.     storeNumLyap = num_lyap;
  577.     num_lyap = 0;
  578.  
  579.     ret = retA; /* if dist below turns out to be small, ret needs a 
  580.             non zero value */
  581.  
  582.  
  583.     /* now the procedure begins */
  584.     while ((dist = distance(ya+zeroth,yb+zeroth,vec_dim)) > IsClose 
  585.         && level == initLevel )
  586.     {
  587.          /* try to refine the pair because they are NOT IsClose */
  588.         /* SET y EQUAL TO THE AVERAGE OF ya AND yb */
  589.         if (BST_switch == 1)
  590.             fraction = (1.0 + random997())/999.;
  591.         else
  592.             fraction = 1./2.;
  593.  
  594.         set_y_fraction(fraction);
  595.         if (ChkKeyStatus() != 0 ) 
  596.         {
  597.             Interrupt();
  598.             if (level < initLevel)
  599.                  continue;/* process terminated via keyboard 
  600.                         Interrupt */
  601.         }
  602.  
  603.         /* Now see where y goes */
  604.         ret = checkTrajectory(level,y);
  605.         scr_rowcol(0,0);
  606.         if (ret == 0)
  607.         {
  608.             erase_line();
  609.             PRINT
  610. "\nret = 0: y didn't go to any attractor; program doesn't know what to do;\n");
  611.             PRINT
  612. "return to single step mode... but we are still in BasinStraddle()        \n");
  613.             cycle = 1;    /* single step mode; */
  614.             Interrupt();
  615.         }
  616.         if (printer == 3) 
  617.         {
  618.             PRINT
  619. "\nab dist=%5.5le, time=%d, dot# =%ld ret=%d (to suppress printing, hit 2)\n",
  620.                             dist,time,dot,ret);
  621.  
  622.         }
  623.         /* now consider case where we fail to refine */
  624.         if (dist < 10.*IsClose && time == 0)
  625.         {    
  626.             PRINT    
  627. "ya[] and yb[] are now in region %d so one has switched basins\n",ret);
  628.             GiveUp();
  629.         }
  630.  
  631.  
  632.  
  633.         set_y_fraction(fraction);/* RECALL y */
  634.  
  635.         /* REFINE ya yb line */
  636.         /* CHOOSE A NEW ya OR yb TO BE THIS MIDPOINT y */
  637.         if (basin[ret] == basinA)
  638.             store(ya,y,lyapzero);/* put y in ya; do not change 
  639.             high coordinates of ya needed for exponents */
  640.         else
  641.         if (ret != 0 && basin[ret] == basinB)/* ret = 0 means the 
  642.                 trajectory did not go to a basin, possibly 
  643.                 because the number of iterates checked = MC 
  644.                 was too low; the routine should stop without 
  645.                 changing yb */
  646.             store(yb,y,lyapzero);
  647.  
  648.     }
  649. /*    if (ret != 0) changed 5/14/87; this is why we gave a default value
  650.         to ret at top of this routine */
  651.  
  652.     OneIterate(ya);
  653.     num_lyap = storeNumLyap;
  654.     OneIterate(yb);
  655. /* Notice that yb[] = y[] now and this is the point that gets plotted and if 
  656.         Lyapunov exponents are being calculated, they are essentially 
  657.         the exponents of yb */
  658. }
  659.  
  660. AccessBS()/* AccessBasinStraddle: this routine 
  661.         takes the pair of vectors ya,yb and refines the pair
  662.         until they are "IsClose"; then it takes a single iterate of
  663.         the refined ya and yb so that the pair advances one iterate;
  664.         it is refined by taking lots of points between ya and yb and 
  665.         as soon as any of them are found to be in basinB, move yb to 
  666.         that point. If none of them are found to be in basinB, choose 
  667.         a point near the midpoint that is known to be a basinA point 
  668.         and move ya to it  */
  669. {
  670.     int initLevel = level;
  671.     int ret, N;
  672.     double dist,distance(),frac,GoldenMean;
  673.     GoldenMean = 2./(1+ sqrt(5.));/* note this is less than 1 */
  674.  
  675.     storeNumLyap = num_lyap;
  676.     num_lyap = 0;
  677.     ret = retA; /* if dist below turns out to be small, ret needs a 
  678.             non zero value */
  679.  
  680.     /* now the procedure begins */
  681.     while ((dist = distance(ya+zeroth,yb+zeroth,vec_dim)) > IsClose 
  682.         && level == initLevel )
  683.     {
  684.          /* try to refine the pair because they are NOT IsClose */
  685.         /* SET y EQUAL TO THE AVERAGE OF ya AND yb */
  686.  
  687.        for (frac = GoldenMean, N = 0; N < divisions && level == initLevel; 
  688.             N++,frac = frac+GoldenMean)
  689.        { 
  690.         while (frac >= 1) 
  691.             frac = frac - 1;
  692.         set_y_fraction(frac);
  693.         /* Now see where y goes */
  694.         ret = checkTrajectory(level,y);
  695.         if (ret == 0)
  696.         {
  697.             scr_rowcol(0,0);
  698.             erase_line();
  699.             PRINT
  700. "\nret = 0: y didn't go to any attractor; program doesn't know what to do;\n");
  701.             PRINT
  702. "return to single step mode... but we are still in BasinStraddle()        \n");
  703.             cycle = 1;    /* single step mode; */
  704.             Interrupt();
  705.         }/* end if */
  706.         if (printer == 3) 
  707.         {
  708.             scr_rowcol(0,0);
  709.             PRINT
  710. "\nab dist=%5.5le, time=%d, dot# =%ld ret=%d (to suppress printing, hit 2)\n",
  711.                             dist,time,dot,ret);
  712.         }/* end if */
  713.         /* now consider case where the intermediate point is close to 
  714.             an attractor (i.e. time == 0 */
  715.         if (dist < 10.*IsClose && time == 0)
  716.         {    
  717.             PRINT    
  718. "ya[] and yb[] are now in region %d so one has switched basins\n",ret);
  719.             GiveUp();
  720.         }/* end if */
  721.  
  722.         /* REFINE ya yb line if basin[ret] == basinB */
  723.         if (basin[ret] == basinB)/* note this includes the case ret = 0
  724.                 so it could cause problems */
  725.         {
  726.             set_y_fraction(frac);/* RECALL y */
  727.             store(yb,y,lyapzero);
  728.             break;/* terminates for() */
  729.         }/* end if */
  730.        }/* end for */
  731.        if (basin[ret] == basinA)
  732.        {
  733.         set_y_fraction(GoldenMean);/* RECALL a  y that gave a 
  734.             ret indicating basinA */
  735.         store(ya,y,lyapzero);
  736.        }
  737.     }/* end while */
  738. /*    if (ret != 0) changed 5/14/87; this is why we gave a default value
  739.         to ret at top of this routine */
  740.  
  741.     OneIterate(ya);
  742.     num_lyap = storeNumLyap;
  743.     OneIterate(yb);
  744. /* Notice that yb[] = y[] now and this is the point that gets plotted and if 
  745.         Lyapunov exponents are being calculated, they are essentially 
  746.         the exponents of yb */
  747. }
  748.  
  749.  
  750. SaddleStraddle()/* we search for a trajectory that stays away from the basins;
  751.          this routine takes the pair of vectors ya,yb and refines the pair
  752.         until they are "IsClose"; then it takes a single iterate of
  753.         the refined ya and yb so that the pair advances one iterate */
  754. {
  755.     int initLevel = level;
  756.     int success;
  757.     char *comment;
  758.     double distance(),RefCallsPerDot,TimePerPlot,dist;
  759.  
  760.     storeNumLyap = num_lyap;
  761.     num_lyap = 0;
  762.     success = 1;/* if ya and yb are close, skip stuff in while and 
  763.             iterate again */
  764.     while ((dist = distance(ya+zeroth,yb+zeroth,vec_dim)) > IsClose 
  765.         && level == initLevel )
  766.     {
  767.             /*try to refine the pair because they are NOT within IsClose */
  768.             success = refineSadStrad();
  769.  
  770.         if (success == 1)
  771.             comment = "success";
  772.         if (success == 2)
  773.             comment = "KEYBOARD TERMINATION";
  774.         if (success == 0)
  775.             comment = "Never exits";
  776.         if (success == -1)
  777.             comment = "refinement failure";
  778.         if (printer >= 2 )
  779.         {
  780.             dist = distance(ya+zeroth,yb+zeroth,vec_dim);
  781.             scr_rowcol(0,0);
  782.             erase_line();
  783.             PRINT
  784. "ab dist=%5.5le, time=%d, dot# =%ld, %s (to suppress printing, hit 2)\n",
  785.                         dist,MaxTime,dot,comment);
  786.             if (dot > 0)
  787.             {
  788.             RefCallsPerDot = totalRefinementCalls/dot;
  789.             TimePerPlot = totalTime/dot;
  790.             erase_line();
  791.             PRINT
  792. "refinement calls per dot = %lf;av. time per plot=%lf;\n"
  793.                 ,RefCallsPerDot,TimePerPlot);
  794.             erase_line();
  795.             PRINT
  796. "totalRefinementCalls=%3.3lf;totalTime=%3.3lf; \n"
  797.             ,totalRefinementCalls,totalTime);
  798.             }
  799.         }
  800.         /* now consider case where we fail to refine */
  801.         if (dist < 10.*IsClose && MaxTime == 0)
  802.         {
  803.             PRINT
  804. "ya[] and yb[] have MaxTime = 0 so process should not proceed\n");
  805.             PRINT
  806. "This suggests points in one ball can actually go to another region\n");
  807.             GiveUp();/* return to single stepping */
  808.         }
  809.     }
  810.  
  811.     if (success == 1)/* this is the good case */
  812.     {
  813.         totalTime = totalTime + time;
  814.         OneIterate(ya);
  815.         num_lyap = storeNumLyap;
  816.         OneIterate(yb);
  817.     }
  818.     num_lyap = storeNumLyap;
  819.   /* Notice that yb[] = y[] now and this is the point that gets plotted and if 
  820.         Lyapunov exponents are being calculated, they are essentially 
  821.         the exponents of yb */
  822. }
  823.  
  824.  
  825.  
  826. int refineSadStrad()
  827.         /* this routine tries to find a refinement, i.e. a new pair
  828.         ya, yb, (one could be unchanged) closer together than the 
  829.         starting pair; if it is successful, the pair is changed and
  830.         1 is returned; 
  831.         it returns 0 if it finds an orbit that never leaves the domain;
  832.         it returns -1 if it cannot find a refinement; notice that all 
  833.         the values of fraction are diadic and so are computed exactly*/
  834. {
  835.     double oldNumLyap;
  836.     double slope, slope0;/* for access. basin strad. */
  837.     double Ytemp[EQUATIONS],fraction=0;
  838.     int initLevel = level;
  839.     int FindTime();
  840.     int numA,numB,monotone,decMonotone,testA;
  841.     int ndivisions; /* a working multiple of divisions */
  842.     int times[1001],N,firstMax;/* firstMax records the first time out of 
  843.         NumSegment+1 times that the value MaxTime is taken on */
  844.  
  845.     totalRefinementCalls++;/* counts the number of refinements */
  846.     for (ndivisions = divisions; ndivisions < 1000
  847.         && level == initLevel ; 
  848.         ndivisions *= 2)/* integer 1000 can't exceed dim times */
  849.     {
  850.         scr_rowcol(6,0);
  851.         MaxTime = 0;
  852.         decMonotone = monotone = 0;
  853.         firstMax = 0;
  854.         for (N = 0; N <= ndivisions && level == initLevel ; N++)
  855.             {
  856. #ifndef MAINFRAME
  857.         if (ChkKeyStatus() != 0 ) 
  858. #endif /* ifndef MAINFRAME */
  859.         {
  860.             set_y_fraction (fraction);
  861.             oldNumLyap = num_lyap;
  862.             Interrupt();/* this may reset the number of Lyapunov 
  863.                 exponens via a #L command; thus we must 
  864.                 check to see if num_lyap has been changed;
  865.                 if so, storeNumLyap is what should have 
  866.                 been set, since num_lyap is supposed to be 
  867.                 0 at this point in the program */
  868.             if (oldNumLyap != num_lyap)
  869.             {
  870.                 storeNumLyap = num_lyap;
  871.                 num_lyap = 0;
  872.             }    
  873.             if (level < initLevel)
  874.                  return (2) ;   /* process terminated via keyboard 
  875.                         interrupt */
  876.         }
  877.         fraction = ((double)N)/ndivisions;/* we want a double 
  878.             precision number */
  879.         FindTime(fraction);/* ret is what checkTra.. returns; 
  880.             "time" is set by FindTime() */
  881.         times[N] = time;
  882.         if (decMonotone == N-1 && times[N-1] >= times[N])
  883.         {
  884.             decMonotone++;/* measures how long the sequence is 
  885.                 non-increasing */
  886.             monotone = decMonotone;/* monotone is supposed to 
  887.                 measure how long times[] is non-decreasing 
  888.                 AFTER the initial period of decrease */
  889.         }
  890.         if (monotone == N-1 && times[N-1] <= times[N])
  891.             monotone++;/* monotone is supposed to 
  892.                 measure how long times[] is non-decreasing 
  893.                 AFTER the initial period of decrease */
  894.  
  895.         
  896.         if (printer >= 3 || gameFlag)
  897.         {
  898.             PRINT    "(%d,%d)  ",N,time);
  899.         }
  900.         if (time > MaxTime)
  901.         {
  902.             firstMax = N;
  903.             MaxTime = time;
  904.         }
  905.         }/* end inner for */
  906.         if (printer >= 2 || gameFlag == YES)
  907.         {
  908.         PRINT 
  909.                 "\n             max point = (%d,%d)\n"
  910.                     ,firstMax,MaxTime);
  911.         if (decMonotone == N-1 && times[0] > times[decMonotone])
  912.             PRINT   
  913.                 "monotone decreasing thru (%d,%d); then    \n",
  914.                 decMonotone,times[decMonotone]);
  915.         else
  916.             PRINT   
  917.                 "constant thru (%d,%d); then      \n",
  918.                     decMonotone,times[decMonotone]);
  919.         PRINT   
  920.             "monotone increasing thru (%d,%d)\n"
  921.                 ,monotone,times[monotone]);
  922.         }
  923. #ifndef MAINFRAME /* cant play interactive game in batch mode */
  924.         if (gameFlag == YES)
  925.         {
  926.             numA = numB = 0;
  927.  
  928.         while (numA == numB)
  929.         {
  930.             PRINT 
  931. "Input two integers >= 0 and <= %d for new ya and yb; NO COMMA between them \n"
  932.                 ,divisions);
  933.             testForInterrupt();
  934.             if (level == initLevel)
  935.             {
  936.                 if (intrptFlag == YES)
  937.                     PRINT 
  938. "Input two integers >= 0 and <= %d for new ya and yb; NO COMMA between them \n"
  939.                         ,divisions);
  940.                 scanf("%d %d",&numA,&numB);
  941.             }
  942.             else
  943.                 break;
  944.             boot = 1;
  945.  
  946.             boot_crt();
  947.         }
  948.         }/* end gameFlag if */
  949. #endif /* ifndef MAINFRAME */
  950.  
  951.     /* ASST: */
  952.         else if (accessFlag == YES)/* decMonotone reports the last in a 
  953.             series of times[] values that do not strictly 
  954.             increase; monotone 
  955.             reports the last subscript in the following series 
  956.             of non decreasing times[] values; hence 
  957.                 monotone > decMonotone;
  958.                                 */
  959.         {
  960.             numA = decMonotone;
  961.         numB = monotone+1;
  962.         if (monotone > ndivisions)
  963.             numB = ndivisions;
  964.         testA = ndivisions/3;
  965.         if (testA < monotone && testA > numA)
  966.         {
  967.             numA = testA;
  968.             if (printer == 3)
  969.                 PRINT 
  970. "choice of numA = ndiv/3 = %d; decMon=%d;\n",numA,decMonotone);
  971.         }
  972.         if (numA == 0 && numB == ndivisions)/* i.e., no change, then 
  973.                        we try a refinement, but it looks bad */
  974.             break;/* break out of for */
  975.         PRINT "Interval refined from (0,%d) to (%d,%d); dot = %ld;\n"
  976.             ,ndivisions,numA,numB,dot);
  977. /*******/
  978.         }
  979.         else /*     SADDLE STRADDLE TRAJECTORY;
  980.         in this case, we find the point in the grid where the time 
  981.         is maximized; then we find the grid points, one on each side 
  982.         which maximize the slope of the times[] function; these grid 
  983.         points are numA and numB, which lets us define ya and yb */
  984.         {
  985.         numA = 0;
  986.         slope0 = 0;
  987.         if (firstMax > 0)
  988.             for (N = 0; N < firstMax; N++)
  989.             {
  990.                 slope = (MaxTime-times[N])/(firstMax-N);
  991.                 if (slope >= slope0)
  992.                 {
  993.                     slope0 = slope;
  994.                     numA = N;
  995.                 }
  996.             }
  997.         else
  998.             numA = -divisions;/* this causes the interval to 
  999.                         expand */        
  1000.         numB = divisions;
  1001.         slope0 = 0;
  1002.         if (firstMax < divisions)
  1003.         {
  1004.             for (N = divisions; N > firstMax; N--)
  1005.             {
  1006.                 slope = (MaxTime-times[N])/(N - firstMax);
  1007.                 if (slope >= slope0)
  1008.                 {
  1009.                     slope0 = slope;
  1010.                     numB = N;
  1011.                 }
  1012.             }
  1013.             if (slope == 0)
  1014.                 numB = 2*divisions;/* this causes the interval 
  1015.                         to expand */        
  1016.         }
  1017.         else
  1018.             numB = 2*divisions;/* this causes the interval to 
  1019.                         expand */        
  1020.         }
  1021.         if (level == initLevel)
  1022.           if (printer >= 2)
  1023.         PRINT "Selected grid points are: %d and %d    \n",numA,numB);
  1024.  
  1025.         /* numA and numB are now set so ya and yb are now changed */
  1026.         fraction = ((double)numB)/ndivisions;    /*we want a double 
  1027.                     precision number */
  1028.         set_y_fraction (fraction);
  1029.         store(Ytemp,y,lyapzero);
  1030.         fraction = ((double)numA)/ndivisions;
  1031.         set_y_fraction (fraction);
  1032.         store(ya,y,lyapzero);
  1033.         store(yb,Ytemp,lyapzero);
  1034.     return (1);/* means refinement is achieved; */ 
  1035.     }    /* end outer for */
  1036.     return (1);/* means refinement is achieved; */ 
  1037. }
  1038.  
  1039.  
  1040.  
  1041.  
  1042.  
  1043.  
  1044. OneIterate(Y)/* iterates the point at Y[] using y */
  1045.     double *Y;
  1046. {
  1047.     store(y,Y,eqn0);
  1048.     (*StraddleIteratee)();
  1049.     store(Y,y,eqn0);
  1050. }
  1051.         
  1052. set_y_fraction(fraction)
  1053.     double fraction;
  1054. {
  1055.     int k;
  1056.     for (k = 0; k < lyapzero; k++) 
  1057.         y[k] = (1-fraction) * ya[k] + fraction * yb[k];
  1058. }
  1059.  
  1060. GiveUp()
  1061. {
  1062.     erase_line();
  1063.     PRINT
  1064.     "Program now enters SINGLE STEP because of this failure.\n");
  1065.     if (storeNumLyap != 0)
  1066.         num_lyap = storeNumLyap;
  1067.  
  1068.     cycle = 1;    /* single step mode */
  1069.     Interrupt();
  1070. }
  1071.  
  1072. int FindTime(fraction)
  1073.     double fraction;
  1074. {
  1075.     int ret;
  1076.  
  1077.     set_y_fraction(fraction);
  1078.     ret = checkTrajectory(level,y);
  1079.     if (ret == 0)        /* == 0 means trajectory did not exit */
  1080.     {
  1081.         set_y_fraction(fraction);
  1082.             /* this will be the value of y that is in y
  1083.              when SaddleStraddle() dies */
  1084.     }
  1085.  
  1086.     return (ret);/*  == 0 if point does not go anywhere */
  1087. }
  1088.  
  1089. triple()
  1090. {
  1091.     double Y[EQUATIONS];
  1092.  
  1093.     set_y_fraction(-1.);
  1094.     store(Y,y,eqn0);
  1095.     set_y_fraction(2.);
  1096.     store(yb,y,lyapzero);
  1097.     store(ya,Y,lyapzero);
  1098. }
  1099.  
  1100. #ifndef MAINFRAME
  1101. testForInterrupt()
  1102.     /*    For detecting non digit when a digit is expected; 
  1103.         any character except \n is put back in the stream "input"; 
  1104.         for anything other than a digit, it returns YES meaning YES 
  1105.         otherwise NO;*/
  1106. {
  1107.     int ch;
  1108.     int ChkKeyStatus();
  1109.     while (ChkKeyStatus() == 0)
  1110.        ;
  1111.     while ((ch = ci())== '\n')/* there may be a return still in 
  1112.         the pipeline from entering the command so this one should 
  1113.         not cause an abort */
  1114.        ;
  1115.     ChkKeyStatus();
  1116.     if (ch >= '0' && ch <= '9')
  1117.     {
  1118.         intrptFlag = NO;
  1119.         ungetc(ch,stdin);
  1120.         return;
  1121.     }
  1122.     else
  1123.     {
  1124.         intrptFlag = YES;
  1125.         PRINT "PAUSE MODE; interrupt is now expected\n");
  1126.         cycle = 1;
  1127.         Interrupt();
  1128.         return;
  1129.     }
  1130. }
  1131. #endif /* ifndef MAINFRAME */
  1132.  
  1133.